Univariate polynomial real root isolation: 
Continued Fractions revisited 

Q I Elias P. Tsigaridas and loannis Z. Emiris 

O ■ Department of Informatics and Telecommunications 

National Kapodistrian University of Athens, HELLAS 
p ' {et , emir is}-(§di .uoa.gr 



(N 



< 



U 



X 



February 1, 2008 



Abstract 



We present algorithmic, complexity and implementation results concerning real root isola- 
tion of integer univariate polynomials using the continued fraction expansion of real algebraic 
C/3 ■ numbers. One motivation is to explain the method's good performance in practice. We im- 

>^, ' prove the previously known bound by a factor of dr, where d is the polynomial degree and 

T bounds the coefficient bitsize, thus matching the current record complexity for real root 
isolation by exact methods. Namely, the complexity bound is OB{d'''T^) using the standard 
^ ' bound on the expected bitsize of the integers in the continued fraction expansion. We show 

^\^ , how to compute the multiplicities within the same complexity and extend the algorithm to 

^^G • non square-free polynomials. Finally, we present an efficient open-source C++ implementation 

'>— f^ ' in the algebraic library SYNAPS, and illustrate its efficiency as compared to other available 

"^ , software. We use polynomials with coefficient bitsize up to 8000 and degree up to 1000. 

5 ! 1 Introduction 

o 



In this paper we deal with real root isolation of univariate integer polynomials, a fundamental 
problem in computer algebra as well as in many applications ranging from computational geometry 
to quantifier elimination. The problem consists in computing intervals with rational endpoints 
}J] ' which contain exactly one real root of the polynomial. We use the continued fraction expansion 

5^ I of real algebraic numbers. Recall that such a number is a real root of an integer polynomial. 

Another motivation is to explain the method's good performance in implementations, albeit 
the higher complexity bounds which was known until now. Indeed, we show that continued 
fractions lead to asymptotic bit complexity bounds that match those recently proven for other 
exact methods, such as Sturm sequences and Descartes' subdivision. 

1.1 Notation 

In what follows Ob means bit complexity and the Os-notation means that we are ignoring loga- 
rithmic factors. For a polynomial A = J2i=i o-i^^ € Z[^], deg {A) denotes its degree. We consider 
square-free polynomials except if explicitly stated otherwise. By C {A) we denote an upper bound 
on the bit size of the coefficients of A (including a bit for the sign). For a € Q, £ (a) > 1 is the 
maximum bit size of the numerator and the denominator. Let M (r) denote the bit complexity 
of multiplying two integers of bit size at most r. Using FFT, M (r) = OsiT'lg'^ t) for a suitable 
constant c. Var{A) denotes the sign variations in the coefficient list of A ignoring zero terms and 
A the separation bound of A, that is the smallest distance between two (complex) roots of A. 
Finally N = max{d, r}. 



1.2 Previous work and our results 

Real root isolation of univariate integer polynomials is a well known problem and various algo- 
rithms exist for it. Moreover, there is a huge bibliography on the problem so we have to mention 
that we only scratch the surface of the existing literature and we encourage the reader to refer to 
the references. 

Exact subdivision based algorithms for real root isolation are based either on Descartes' rule of 
sign or on Sturm sequences. Roughly speaking, the idea behind both approaches is to subdivide a 
given interval that initially contains all the real roots until it is certified that none or one root is 
contained. Descartes' approach achieves this by repeatedly transforming the original polynomial 
and counting the sign variations in the coefhcients' list, while Sturm's approach computes a signed 
polynomial remainder sequence and evaluates it over the endpoints of the interval of interest. 
Quite recently it was proven (cf ^1 E| and references therein) that both approaches, the one 
based on Descartes' rule of sign (where the polynomials are represented either in the monomial 
or in the Bernstein basis) and the one based on Sturm sequences achieve the same bit complexity 
bound, namely Osid'^T^) or Ob{N^)- Moreover using Sturm sequences in a pre-processing and a 
post-processing step |21| the bound holds for the non square-free case and the multiplicities of the 
roots can also be computed. 

The continued fraction algorithm (from now on called CF) differs from the subdivision based 
algorithms in that instead of bisecting a given initial interval it computes the continued fraction 
expansions of the real roots of the polynomial. The first formulation of the algorithm is due 
to Vincent 001) see also 2 for historical references, based on his theorem (Th. 0] without the 
terminating condition) where it was stated that repeated transformations of the polynomial will 
eventually yield a polynomial with zero (or one) sign variation, thus Descartes' rule implies the 
transformed polynomial has zero (resp. one) real root in (0,oo). If one sign variation is attained 
then the inverse transformation can be applied in order to compute an isolating interval for the real 
root that corresponds to the original polynomial and moreover the q's appear in the transformation 
correspond to the partial quotients of the continued fraction expansion of the real root. However 
Vincent's algorithm is exponential ll^l . He computed the q's in the transformation of Th. 0]by 
repeated shift operations of the form X i— > X + 1, thus if one of the q's (or even the sum of all) 
is of magnitude, say, 2"^ then an exponential number of steps must be performed. 

Uspensky ;37) extended Vincent's theorem by computing an upper bound on the number of 
transformations so as to isolate the real roots, but failed to deal with its exponential behavior. See 
also '12' '32^ where the problem of approximating a real algebraic number is also considered. Using 
Vincent's theorem, Collins and Akritas J_3 derived a polynomial subdivision-based algorithm using 
Descartes' rule of sign. Akritas [HI Q dealt with the exponential behavior of the CF algorithm, 
by computing the q's in the transformations as positive lower bounds of the positive real roots, 
via Cauchy's bound (for details, see sec. O. He achieved a complexity of Osid^T^) or Ob{N^), 
without using fast Taylor shifts ^1^. However, it is not clear how this approach accounts for the 
increased coefficient size in the transformed polynomial after applying X i— > b + X. Another issue 
is to bound the size of the Ci . Refer to Eq. Q which indicates that the magnitude of the partial 
quotients is unbounded. CF is the standard real root isolation function in Mathematica 4 and 
for some experiments against subdivision-based algorithms, also in Mathematica, the reader 
may refer to 0. 

Another class of univariate solvers are numerical solvers, e.g. PHI IHl IH] that compute an 
approximation of all the roots (real and complex) of a polynomial up to a desired accuracy. 

The contributions of this paper are the following: First, we improve the bound of the number 
of steps (transformations) that the algorithm performs. This is basically achieved through Th.|S| 
Second, we bound the bitsize of the partial quotients and thus the growth of the transformed 
polynomials which appear during the algorithm. For this we use the theory of the continued 
fraction expansion of real numbers and a standard average case analysis. We revisit the proof of 
ISm so as to improve the overall bit complexity bound of the algorithm to Ob{N^), thus matching 
the current record complexity for real root isolation. The extension to the non square-free case 
uses the techniques from |18| . Third, we present our efficient open-source C++ implementation 



and illustrate it on various data sets, including polynomials of degree up to 1000 and coefficients 
of 8000 bits. Our software seems faster than the root-isolation implementations that we tested, 
including RS. We also tested a numeric solver, namely ABERTH, which has comparable efficiency 
and, on many instances, is slower. We believe that our software contributes towards reducing the 
gap between rational and numeric computation, the latter being usually perceived as faster. 

The rest of the paper is structured as follows. The next section sketches the theory behind 
continued fractions. Sec. O presents the CF algorithm and Sec. 0] its analysis. We conclude with 
experiments using our implementation, along with comparisons against other available software 
for univariate equation solving. 

2 Continued fractions 

We present a short introduction to continued fractions, following |38| which although is far from 
complete suffices for our purposes. The reader may refer to e.g [51 I42L [TUl I38| . In general a simple 
(regular) continued fraction is a (possibly infinite) expression of the form 

Co H = [co,Ci,C2,...] 

Cl 



C2 + . . . 

where the numbers Ci are called partial quotients, Ci e Z and Cj > 1 for i > 0. Notice that cq may 
have any sign, however in our real root isolation algorithm cq > 0. By considering the recurrent 
relations 



P-l = 1, Po^ Co, Pn+1 = C„+i P„ + Pn-1 




Q-1 = 0, Qo = 1, Qn+l = Cn+lQn + Qn-l 




it can be shown by induction that i?„ = £^ = [cq, ci, . . . , c„], for n = 0, 1, 2, . 


. and moreover 


that 

PnQn+l-Pn+lQn = (-!)"+! 




PnQn+2-Pn+2Qn = (-l)"+lc„+2 





If 7 = [co, Cl,...] then 7 = co + q^ - q^ + ■ ■ ■ = cq + ^^^^^ qSiQ,. ^^^ ^™*^^ ^^^^ 
is a series of decreasing alternating terms it converges to some real number 7. A finite section 
Rn = Q^ = [co, Cl, . . . , c„] is called the n— th convergent (or approximant) of 7 and the tails 
7n+i = [cn+1, c„+2, • ■ ■ ] are known as its complete quotients. That is 7 = [co, ci, . . . , c„, 7„+i] for 
n = 0, 1, 2, . . . . There is a one to one correspondence between the real numbers and the continued 
fractions, where evidently the finite continued fractions correspond to rational numbers. 

It is known that Qn > Fn+i and that Fn+i <</)"< -^^+2, where Fn is the n— th Fibonacci num- 
ber and (f> = ^2^ ^^ ^^^ golden ratio. Continued fractions are the best (for a given denominator 
size), approximations, i.e 



1 

< 



]n{Qn+l + Qn) 



Ik 

^71 



< — < T^ < 

1n+l 



QnQr. -^ ^ n2 



Let 7 = [cojCi,...] be the continued fraction expansion of a real number. The Gauss-Kuzmin 
distribution [101 12] states that for almost all real numbers 7 (meaning that the set of exceptions 
has measure zero) the probability for a positive integer 5 to appear as an element in the continued 
fraction expansion of 7 is 

Pro5[c, = <5]=lg|^ti|, ,>o (1) 

The Gauss-Kuzmin law induces that we can not bound the mean value of the partial quotients 
or in other words that the expected value (arithmetic mean) of the partial quotients is diverging, 
E[ci\ — X^i^i ^ P^oh[ci = 5] = 00, z > 0. Surprisingly enough the geometric (and the harmonic) 



Algorithm 1: CF(A,M) 



Input: AeZ[X],M(X) = ;^f±^, k,l,m,neZ 
Output: A list of isolating intervals 

1 if A{0) = then 

2 OUTPUT Iiiterval( M(0),M(0)) ; 

3 .4^.4(A)/A; 

4 CF(A,M); 

5 V ^ Var(A); 

6 if F = then return ; 

7 if 1/ = 1 then 

8 OUTPUT Interval ( Af (0), Af (oo)); 

9 RETURN ; 

10 6 ^ PLB(yl) // PLB = PositiveLowerBound ; 

11 if 6 > 1 then A^ A{b + A), A/ ^ Af (6 + A) ; 

12 Ai ^ A{1 + A), Ml ^ Af (1 + A) ; 

13 CF(Ai,Afi) // Looking for real roots in (l,+oo) 

14 A2 ^ A(^), Af2 ^ Af(^) ; 

15 CF{A2,M2) II Looking for real roots in (0,1) ; 

16 RETURN : 



mean is not only asymptotically bounded, but is bounded by a constant. For the geometric mean 
this is the famous Khintchine's constant ^23j, i.e. 



lim 



\ 



J|ci=/C = 2.685452001. 



i=\ 



which is an irrational number. The reader may refer to [H] for a comprehensive treatment of 
Khintchine 's means. The expected value of the bitsize of the partial quotients is a constant for 
almost all real numbers, when n —^ oo or n sufficiently big |23ll31j . Following closely \61\ . we have: 
£:[lnc,] = ^ELilnc, = In/C = 0.98785..., as n -> oo, Vi > 0. Let £(c,) = b„ then 

m] = Oil) (2) 

A real number has an (eventually) periodic continued fraction expansion if and only if it is a root 
of an irreducible quadratic polynomial. "There is no reason to believe that the continued fraction 
expansions of non-quadratic algebraic irrationals generally do anything other than faithfully follow 
Khintchine's law" ^J, and also various experimental results |1U[I31|I5^ suggest so. For the largest 
digit that can appear in the partial quotients of a rational number the reader may refer to [221 ■ 

3 The CF algorithm 

Theorem 1 (Descartes' rule of sign) The number R of real roots ofA{X) in (0,oo) is bounded 
by Var{A) and we have R = Var{A) mod 2. 

Remark 2 In general Descartes' rule of sign obtains an overestimation of the number of the 
positive real roots. However if we know that A is hyperbolic, i.e has only real roots or when the 
number of sign variations is or 1 then it counts exactly. 

Theorem 3 (Budan) }26l \^ Let a polynomial A, such that deg(A) — d and let a < b, where 
a,b gM.. Let Aa, resp. Ab, be the polynomial produced after we apply the map A i— > A + a, resp. 



X 1-^ X + b, to A. Then the fallowings hold: (i) Var(Aa) > Var{Ab), (ii) #{7 G (a, 6)|A(7) = 
0} < Var{Aa) - VariAt) and (iii) #{7 € (a, b)\A{-i) = 0} = Var{Aa) - Var{Ab) mod 2. 

The CF algorithm depends on the following theorem, which dates back to Vincent's theorem 
in 1836 HOI- The inverse of Th.Hthat proves the termination of CF can be found in pi lTll^ . 
It is a very interesting question whether the one and two circle theorems (c.f ^Sl and references 
therein), employed in the analysis of the subdivision-based real- root isolation algorithm |13| . can 
also be applied and possibly improve the complexity of the CF algorithm. 

Theorem 4 |^^ Let A € Z[X], with deg(yl) — d and let A > be the separation bound. Let 
n be the smallest index such that 

F„_iA>2 and F„_iF„A > IH 

(d 

where Fn is the n-th Fibonnaci number and e^ = (1+^) ^i^^— 1. Then the map X i— > [cq, Ci, . . . , c„, X], 
where Cg, Ci, . . . , c„ is an arbitrary sequence of positive integers, transforms A{X) to An{X), which 
has no more than one sign variation. 

Remark 5 Since -^ < £d < -§1 ',14^ we conclude that ^ -I- 1 < 2d^ for d> 2. Thus, if d> 2 we 
can replace the two conditions of Th. 0] by Fn-iA > 2d^ , since Fn > Fn-i > 1 and Fn-iFnA > 
F„_iA > 2^2 > 2. 

Th. 0] can be used to isolate the positive real roots of a square-free polynomial A. In order to 
isolate the negative roots we perform the transformation X 1-^ —X, so in what follows we will 
consider only the positive real roots of A. 

Vincent's variant of the CF algorithm goes as follows: A polynomial A is transformed to Ai by 
the transformation X 1-^ 1+X and if Var{Ai) = or Var{Ai) = 1 then A has 0, resp. 1, real root 
greater than 1 (Th.0. IiVar{Ai) < yar(A) then (possibly) there are real roots of A in (0, 1), due 
to Budan's theorem (Th.|21l. A2 is produced by applying the transformation X 1— > 1/(1 -I- X) to 
A, if Var{A2) = or Var{A2) = 1 then A has 0, resp. 1, real root less than 1 (Th.^. Uspensky's 
pTZJ variant of the algorithm (see also [22]) at every step produces both polynomials Ai and A2. 
probably, as Akritas states l2j, because he was unaware of Budan's theorem (Th. |3J|. In both 
variants, if the transformed polynomial has more than one sign variations, we repeat the process. 

We may consider the process of the algorithm as an infinite binary tree in which the root 
corresponds to the original polynomial A. The branch from a node to a right child corresponds 
to the map X h^ AT + 1, while to the left child to the map X h^ j+x- Notice that a sequence of 
c transformations A 1— > 1 -I- A followed by one of the type A 1-^ 1/(1 + A) is equivalent to two 
transformations, one of the type A 1— > c -I- 1/X followed by A i-^ 1 + A. Thus Vincent's algorithm 
(and Uspensky's) results to a sequence of transformations like the one described in Th. 01 and 
so the leaves of the binary tree that we considered hold (transformed) polynomials that have no 
more than one sign variations, if Th. ^ holds. Akritas E] replaced a series of A 1— > A -I- 1 
transformations by A i-^ A + 6, where b is the positive lower bound (PLB) on the positive roots of 
the tested polynomial. This was computed by Cauchy's bound O |^ 021 ■ This way, the number 
of steps is polynomial and the complexity is in Osid^T^)- However, it is not clear whether or how 
the analysis takes into account that the coefficient bitsize increases after a shift. Another issue is 
to bound the size of the q. 

For these polynomials that have one sign variation we still have to find the interval where the 
real root of the initial polynomial A lies. Consider a polynomial A„ that corresponds to a leaf 
of the binary tree that has one sign variation. Notice that A„ is produced after a transformation 
as in Th. ^ using positive integers cq, ci, . . . , c„. This transformation can be written in a more 
compact form using the convergents 

QnA -I- Qn~l 



where „" ^ and jf- are consecutive convergents of the continued fraction [cq, ci, . . . , c„]. Notice 
that |(31l is a Mobius transformation, see 01121 for more details. Since An has one sign variation 
it has one and only one real root in (0, oo), so in order to obtain the isolating interval for the 
corresponding real root of A we evaluate the right part of Eq. |(2Jl once over and once over cxd. 
The (unordered) endpoints of the isolating interval are jf—^ and ^ . 

The pseudo-code of the CF algorithm is presented in Alg.Q] Notice that the Interval function 
orders the endpoints of the computed isolating interval and that PLB(A) computes a lower bound 
on the positive roots of A. The initial input of the algorithm is a polynomial A{X) and the 
trivial transformation M{X) = X. We need the functional M in order to keep track of the 
transformations that we perform so that to derive the isolating intervals. Notice that Lines ^ and 
^are to be executed only when Var{Ai) < Var{A2), but in order to simplify the analysis we omit 
this, since it only doubles the complexity. 

4 The complexity of the CF algorithm 

The complexity of the algorithm depends on the number of transformations and the cost of each. 
However special care should be taken since after each transformation the bit size of the polynomial 
increases. 

Let disc(A) be the discriminant and lead (A) the leading coefficient of A. Mahler's measure 
of a polynomial A is A4{A) = | lead (A) | ni=i rnaxjl, |7i|}, where 7^ are all the (complex) roots 
of A |7||32Enil2Z|- Moreover M{ A) < T^d+ 1. We prove the following theorem, which is based 
on a theorem by Mignotte |5H], thus extending [T^ITt] . 

Theorem 6 Let A E 1\X\, with deg{A) — d and £ (A) = t. Let ft be any set of k couples of 
indices {i,j) such that 1 < i < j < d and let the non-zero (complex) roots of A be < \"/i\ < I72I < 
• • • < \jd\- Then 

2''MiA)''> n l7» -lj\> 2'=-^^^ MiAf-""-^ y/disc{A) 
(jj)Gn 

Proof. Let il be the multiset H. — {j\{i,j) G ft} and \n\ ~ k. We use the inequality |a — 5| < 
2max{|a|,|6|}(*)a,6e C and the fact EHIEZI that for any root of A, j^;^ < |7,;| < MiA). 
In order to prove the left inequality 

n l7. - 7,1 < 2'^ n l^^l ^ 2'= ma^ |7,f < 2>^M{A)\ 

Recall 131111123 that disc(A) = lead(A)^'^"^ JJ^^^ (7^ -7^)^. For the right inequality we 
consider the absolute value of the discriminant of A: 

|disc(A)| = \leB.diA)\'^~'U^<,h~lj? 

= \le^.<iiA)\''^-'Ui.,,)en\l^~l,\'Ui.,Mn\l^-^^\' ^ 
v/|disc(A)| = I lead (A) l^-^ Ui^,Mi, h^ " lj\ Ui^,Mn h^ " ^j\ 

We consider the product Yiu j)dn l7« ~ 7il ^^'^ ^*^ apply ^ ^ ' — k times inequality (*), thus 



n(,,)^al7.-7,l < 2^^-'hini2\'---h,\'-'iU,enhA)-' .4^ 

< 2^^^-''M{A)d-^\lea.d{A)\^-''MiA)'' 

where we used the inequality |7i|°|72|^ • • • \ld\'^^^ < \M{A)/ le&d{A) \'^-^, and the fact ^ that, 
since Vi, \ji\ > M{A)^^, we have n,en l7il - l7il'^ - ■M{A)^''. Thus wc conclude that 



Uu.mn\l^-lj\ > 2^-^^Al(A)i-'^-Vldisc(A)| 



n 

A similar theorem but with more strict hypotheses on the roots first appeared in "IS" and the 
conditions were generalized in J17j : namely in order for the bound 15^, 17 to hold the sets of 
indices i and j should be rearranged such that they form an acyclic graph where each node has 
out-degree at most one. The bound of Th. Elhas an additional factor of 2'^ wrt ^|E], which 
plays no role when d = 0{t) or when notation with N is used. Moreover we loosen the hypotheses 
of the theorem and thus all the proofs concerning the number of steps of the subdivision-based 
solvers |17II21| are dramatically simplified. Possibly a more involved proof of Th.Elmay eliminate 
this factor -'^. 

Remark 7 There are two simple however crucial observations about T/i.^j When the transformed 
polynomial has one sign variation, then the interval with endpoints "''" — [cq,Ci, . . . ,c„_i] and 

Wn — 1 

^ = [cqjCi, . . . , c„] (possibly unordered) isolates a positive real root of A, say 7^. Then, in order 
for Th. P] to hold, it suffices to consider, instead of the separation bound A, the quantity [7^ — 7ci |, 
where ^^ is the (complex) root of A closest to 7^. When the transformed polynomial has no sign 
variation and [cq,Ci,. . . ,c„] is the continued fraction expansion of the (positive) real part of a 
complex root of A, say "fi, then again it suffices to replace A by \ji — 7cJ. 

Theorem 8 The CF algorithm performs at most 0{d^ + dr) steps. 

Proof. Let < I71I < I72I < • • • < |7fe|, fc < d be the (complex) roots of A with positive real part 
and let ^a denote the root of A that is closest to 7^. 

We consider the binary tree T generated during the execution of the CF algorithm. The number 
of steps of the CF algorithm corresponds to the number of nodes in T, which we denote by #(r). 
We use some arguments and the notation from |17j in order to prune the tree. 

With each node w of T we associate a Mobius transformation M^ : X 1-^ ^^~^n ' ^ polynomial 
Ay and implicitly an interval ly whose unordered endpoints can be found if we evaluate M„ on 
and on 00. Recall that A^ is produced after M^ is applied to A. The root of T is associated with 
A, M{X) = X (i.e fc = n=l,/ = m = 0) and implicitly with the interval (0, 00). 

Let a leaf w of T be of type-i if its interval /„ contains i > real roots. Since the algorithm 
terminates the leaves are of type-0 or type-1. We will prune certain leaves of T so as to obtain 
a certain subtree T' where it is easy to count the number of nodes. We remove every leaf that 
has a sibling that is not a leaf. Now we consider the leaves that have a sibling that is also a leaf. 
If both leaves are of type-1, we arbitrary prune one of them. If one of them is of type-1 then we 
prune the other. If both leaves are of type-0, this means that the polynomial on the parent node 
has at least two sign variations and thus that we are trying to isolate the (positive) real part of 
some complex root. We keep the leaf that contains the (positive) real part of this root. And so 
#(T)<2#(T'). 

Now we consider the leaves of T' . All are of type-0 or type-1. In both cases they hold the 
positive real part of a root of A, the associated interval is |/i,| > [7^ — ^a \ (Rem.Ej) and the number 
of nodes from a leaf to the root is n^, which is such that the conditions of Rem. are satisfied. 
Since Ui is the smallest index such that the condition of Rem. hold, if we reduce Ui by one then 
the inequality does not hold. Thus 

Fn^-2\l^ - 7cJ < 2^2 ^ 0"'-3|7, - 7,J < 2^2 ^ n, < 4 + 21gd- lg|7, - 7,J 

We sum over all n^ to bound the nodes of T', thus 

#(r') < ^ n» < 2fc(2 + Igd) - ^ log |7. - 7cJ < 2fc(2 -I- Igd) - logf[ I7, - 7,, | (5) 

i— 1 z— 1 i— 1 

In order to apply Th. we should rearrange Jli=i 17* ~ lci\ ^^ that the requirements on the 
indices of roots are fulfilled. This can not be achieved when symmetric products occur and thus the 

^Personal communication with M. Mignotte 



worst case is when the product consists only of symmetric products i.e YiiLi lilj ~ 7c )(7c — 7j)l- 
Thus we consider the square of the inequahty of Th.|H| taking | instead of k and disc(yl) > 1 
(since A is square- free), thus 

ntil7.-7c.| > {2^^"^ MiA^-'-^y ^g^ 

-logntil7.-7cJ < d^~d-k+{2d + k-2)\gM{A) 

Eq. (0 becomes #(T') < 2fc(2 + Igd) + d^ - d- k + {2d + k - 2)\gM{A). However for 
Mahler's measure it is known that M{A) < 2''^/d+ 1 ^ IgA^(A) < r + Igd, for d > 2, thus 
#(r') < 2k{2 + lgd) + d^ -- d - k + {2d + k - 2){t + \gd). Since #(T) < 2#(T') and k < d, we 
conclude that # (T) = O (d^ + rf r + d Ig d) . D 



4.1 Real root isolation 

To complete the analysis of the CF algorithm we have to compute the cost of every step that the 
algorithm performs. In the worst case every step consists of a computation of a positive lower 
bound b (Linenj and three transformations, X h^ b+X , X h^ 1 + X and X i— > yqi^ (Linesnjnand 
^in Alg.^l. Recall, that inversion can be performed in 0{d). Thus the complexity is dominated 
by the cost of the shift operation (Linen in Alg.QJ if a small number of calls to PLB is needed in 
order to compute a partial quotient. We will justify this in Sec. 14.21 In order to compute this cost 
a bound on C{ck) — bk,0 < k < rui is needed, see Eq. ^. 

For the analysis of the CF algorithm we will need the following: 

Theorem 9 (Fast Taylor shift) gT)/ Let AeZ[X], with deg{A) ^ d and C {A) = r and let a e 

Z, such thatC (a) = a. Then the cost of computing B — A{a+X) G Z[X] isOsiM {d^lgd + d^a + dr)). 

Moreover C{B) =0{t + da). 

Initially A has degree d and bitsize t. Evidently the degree does not change after a shift 
operation. Each shift operation by a number of bitsize bh increases the bit size of the polynomial 
by an additive factor dbh^ in the worst case (Th. |5J|. At the /i— th step of the algorithm the 
polynomial has bit size 0{t + dX]i=i ^i) ^^'^ ^^ perform a shift operation by a number of bit 
size &;i+i. Th. 121 states that this can be done in Ob ( M (d^\gd + d^b^+i + d{T + dX]i=i ^«) ) ) o^ 
Ob (m (d2 \gd + dT + d^ E'=i' h) ) . 

Now we have to bound X]i=i ^i- F'or this we use Eq. (0), which bounds E[bi]. By linearity 
of expectation it follows that E[J2'l;^^ bi\ = 0{h) Since h < #(T) = 0{d'^ + dr) (Th. ^, the 
(expected) worst case cost of step h is Ob{M {d^ \gd + dr + d'^{d'^ + dr))) or OB{d^{d'^ + dr)). 
Finally, multiplying by the number of steps, 4t^{T), we conclude that the overall complexity is 
dB{d'^ + d^T + d'^T'^), OTOBid^T^) iid = 0{T), ot Ob {N^), where N = max {d,T}. 

Now let us isolate, and compute the multiplicities, of the real roots oi Ain € Z[X], which is not 
necessarily square- free, with deg(Ai„) = d and C {Ain) — t. We use the technique from [21] and 
compute the square- free part A of Ain using Sturm-Habicht sequences in OB{d^T). The bit size 
of A is £ {A) = 0{d + T). Using the CF algorithm we isolate the positive real root of A and then, 
by applying the map X i~> ~X, we isolate the negative real roots. Finally, using the square-free 
factorization oi Ain, which can be computed in ©^(d'^T), it is possible to find the multiplicities in 
dB{d^T). 

The previous discussion leads to the following theorem. 

Theorem 10 Let A £ Z[X] (not necessarily square-free) such that deg,{A) = d > 2 and C {A) = r. 
We can isolate the real roots of A and compute their multiplicities in expected time OB{d^ -l-d^r^), 
or dB{N^), where N = max {d,T}. 

Remark 11 The hypothesis d> 2 may be replaced by d> A, since real solving of polynomials of 
degree up to 4 can be performed in 0(1) or Ob{t) }lf^ - 



4.2 Rational roots and PLB (Positive Lower Bound) realization 

This section studies a way to compute a lower bound on the positive roots and presents its efRciency 
and accuracy. It seems that this is the standard approach in CF algorithms, though it is seldom 
discussed. 

Let us consider the special case of rational roots. Their continued fraction expansion does 
not follow Khintchine's law, even though there are results |22] on the largest digit in such an 
expansion. However, recall that if ^ is a root of A then p divides Oq and q divides a^j, thus in 
the worst case £ {p/q) = 0{t) and so the rational roots are isolated fast. Treating them as real 
algebraic numbers leads to an overestimation of the number of iterations. 

There is one exception to this good behavior of rational roots, namely when they are very large, 
well separated, and we are interested in practical complexity ^. This is due to the fact that PLB 
must be applied many times. In 31 , the authors performed a small number of Newton iterations 
in order to have a good approximation of a partial quotient. In 01|3|, this problem was solved 
by performing a homothetic transformation, X i— > bX, where b is the computed bound whenever 
b > 16. We follow the latter approach so, after Line Q] in Alg. ^ if 6 = PLB(yl) > 16, we apply 
X K^ bX. The cost of applying X i— > bX is in the worst case the same as doing a shift. 

Now, let us see how PLB(A) is obtained in general. It is computed as the inverse of an upper 
bound on the roots of X'^A{^). In general PLB(A) is applied more than once in order to compute 
some Ci. However this number is very small jSIH- Eq- iQ) implies that the probability that a 
partial quotient is < 10 is ~ 0.87, thus in general the partial quotients are of small magnitude. 
Moreover it is known |10[I38| that 7 ^ §^ < - — ^-qt, which means that an extremely big value of 
a partial quotient indicates that the previous approximation of the algebraic number was extremely 
good, thus it will need a small number of steps to isolate it. 

This is how we implement PLB: Let [/ = {j G N | j < d A a^ < 0} and assume a^ > 0. Now let 
C = adx'^ + SiG(7 ^i^^ ^ Z[j\r]. Then C has exactly one nonnegative root, say t, which is an upper 
bound on the roots of A. We set PLB(A) = 2maxa.gc/ |.2i|i/i^ which is nearly optimal |21], or to 
be more specific t < PLB(yl) < 2t. Actually this bound "[...] is to be recommended among all" |39| . 
since it provides a very good approximation and is easily implementable. In our implementation 
we compute PLB only as powers of 2 so that we can take advantage of fast operations as in |SS| ■ 
Notice that the cost of computing PLB is small, namely 0{d) 0. 

We can improve the computed bound by applying a small number of a modified Newton's 
iteration [^S] |^ to C, that is guaranteed to converge rapidly. Notice that PLB is not a general 
bound on the roots, like e.g. |42[ I26L 1^ 1271 134| . but a bound on the positive roots only, see |24II36| . 
Moreover, when the number of negative coefficients is even then a bound due to Stefanescu |36| 
can be used which is much better. 

5 Implementation and experiments 

We have implemented the CF algorithm in SYNAPS ^ ,22,, which is a C++ library for symbolic- 
numeric computations that provides data-structures, classes and operations for univariate and 
multivariate polynomials, vector, matrices, ... Our code will be included in the next public release 
of SYNAPS. The implementation is based on the integer arithmetic of GMP'^ (v. 4.1.4) and uses 
only transformations of the form X i-^ 2^X and X i-^ X + 1 to benefit from the fast implemen- 
tations that are available in GMP. However, our implementation follows the generic programming 
paradigm, thus any library that provides arbitrary precision integer arithmetic can be used instead 
of GMP. 

We restrict ourselves to square-free polynomials of degree € {100, 200, . . . , 1000}. Following 
|33|. the first class of experiments concerns well-known ill-conditioned polynomials namely: La- 
guerre (L), first (CI) and second (C2) kind Chebyshev, and Wilkinson (W) polynomials. We also 
consider Mignotte (Ml) polynomials X"^ — 2(101X — 1)^, that have 4 real roots but two of them 

^www-sop. inria.fr/galaad/logiciols/synaps/ 
3www.swox.coni/gmp/ 



very close together, and a product (X'' - 2(101X - l)^) (X'^ - 2((101 + j^)X - l)^) of two such 
polynomials (M2) that has 8 real roots. Finally, we consider polynomials with random coefhcients 
(Rl), and monic polynomials with random coefficients (R2) in the range [-1000, 1000], produced 
by MAPLE, using 101 as a seed for the pseudo-random number generator. 

We performed experiments against RS "*, which seems to be one of the fastest available software 
for exact real root isolation. It implements a subdivision-based algorithm using Descartes' rule of 
sign with several optimizations and symbolic- numeric techniques [IM^\ . Note that we had to use RS 
through its MAPLE interface. Timings were reported by its internal function rs_time(). 

We also test ABERTH |Hlini> which a numerical solver with unknown (bit) complexity but very 
efficient in practice, available through SYNAPS. In particular, it uses multi-precision floats and 
provides a floating-point approximation of all complex roots. Unfortunately, we were not always 
able to tune its behavior in order to produce the correct number of real roots in all the cases, i.e. 
to specify the input and the output precision. 

In SYNAPS, there are several univariate solvers, based on Sturm sequences, Descartes' rule of 
sign, Bernstein basis, etc (see [18 for details and experimental results). CF is clearly faster than 
all these solvers, therefore we do not report on these experiments. In particular, the large inputs 
used here are not tractable by the Sturm-sequence solver in SYNAPS, and this is also the case for 
another implementation of the Sturm-sequence solver in CORE ^ . 

So, in TableQ] we report experiments with CF, RS, ABERTH, where the timings are in seconds. 
The asterisk (*) denotes that the computation did not finish after 12000s and the question-mark 
(?) that we were not able to tune the ABERTH solver. The experiments were performed on a 
2.6 GHz Pentium with 1 GB RAM, and our code was compiled using g++ 3.3 with options -03 
-DNDEBUG. 

For (Ml) and (M2), there are rational numbers with a very simple continued fraction expansion 
that isolate the real roots which are close. These experiments are extremely hard for RS. On (Ml), 
ABERTH is the fastest and correctly computes all real roots, but on (M2), which has 4 real roots 
close together, it is slower than CF. CF is advantageous on (W) since, as soon as a real root is 
found, transformations of the form X i-^ X + 1 rapidly produce the other real roots. We were 
not able to tune ABERTH on (W). For (L), (CI) and (C2), CF is clearly faster than RS, while we 
were not able to appropriately tune ABERTH to produce the correct number of real roots. The 
polynomials in (Rl) and (R2) have few and well separated real roots, thus the semi-numerical 
techniques in RS are very effective. To be more specific, RS isolates all roots using only 63 bits of 
accuracy (this information was extracted using the function rs_verbose( 1)). ABERTH is even 
faster on these experiments. However, even in this case, CF is faster than RS; it is a little slower 
than ABERTH (see Table ^ . 

We finally tested a univariate polynomial that appears in the Voronoi diagram of ellipses [201 ■ 
The polynomial has degree 184, coefhcient bitsize 903, and 8 real roots. CF solves it in 0.12s, RS 
in 0.3s and ABERTH in 1.7s. 

There are ways to improve our solver. First, instead of exact integer arithmetic we may 
use semi- numerical techniques like those in RS j33| . These techniques may be based on interval 
arithmetic. 
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